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ABSTRACT 



Context. When a neutron star's rotation slows down, its internal density increases, causing deviations from beta equilibrium that 
induce reactions, heating the stellar interior. This mechanism, named rotochemical heating, has previously been studied for non- 
superfluid neutron stars. However, the likely presence of superfluid nucleons will affect the thermal evolution of the star by suppress- 
ing the specific heat and the usual neutrino-emitting reactions, while at the same time opening new Cooper pairing reactions. 
Aims. We describe the thermal effects of Cooper pairing with spatially uniform and isotropic energy gaps of neutrons A„ and protons 
A p , on the rotochemical heating in millisecond pulsars (MSPs) when only modified Urea reactions are allowed. In this way, we are 
able to determine the amplitude of the superfluid energy gaps for the neutron and protons needed to produce different thermal evolu- 
tion of MSPs. 

Methods. We integrate numerically, and analytically in some approximate cases, the neutrino reactions for the modified Urea pro- 
cesses with superfluid nucleons to include them in the numerical simulation of rotochemical heating. 

Results. We find that the chemical imbalances in the star grow up to the threshold value A, hr = min^A,, + 3A p , 3A„ + A p ), which 
is higher than the quasi-steady state achieved in the absence of superfluidity. Therefore, the superfluid MSPs will take longer to 
reach the quasi-steady state than their nonsuperfiuid counterparts, and they will have a higher a luminosity in this state, given by 
Ly' qs = (1 - 4) x 10 32 (A,i, r /MeV) (P-w/Pms) er 8 s wnere P-20 is the period derivative in units of 10~ 20 and f ms is the period in 
milliseconds. We are able to explain the UV emission of the PSR J0437-4715 for 0.05[MeV] < A thr < 0.45[MeV]. These results are 
valid if the energy gaps are uniform and isotropic. 

Key words, stars: neutron — dense matter — relativity — stars: rotation — pulsars: general — pulsars: individual (PSR J0437-4715) 



1. Introduction 

The observation of thermal emission from the surface of a 
neutron star (NS) has the potential to provide constraints on 
its inner structure. In the existing literature, several detailed 
cooling calculations have been compared to the few estimates 
available for the surface temperatures of neutron stars (see 
I Yakovlev & Pethickl I2004I for a review and references). These 
calculations are based on passive cooling, at first neutrino- 
dominated, and later driven by photon emission at the age of 
~ 10 5 yr. 

Several mechanisms can keep NSs hot beyond the standard 
cooling timescale ~ 10 7 yr. One of them is rotochemical heat- 
ing, which has previously been studied for neutron stars of non- 
superfluid matter. It was first proposed in iReiseneggerf (1995) 
and then developed in lFernandez & Reisene gger (2005) by con- 
sidering the internal structure via realistic EOSs, in the frame- 
work of general relativity. It works as follows. The reduction in 
the centrifugal force makes the NS contract. This perturbs each 
fluid element, raising the local pressure and causing deviations 
from beta equilibrium. The system eventually reaches a new 
quasi-equilibrium configuration where the rate at which spin- 
down modifies the equilibrium concentrations is the same as that 
at which neuttino reactions restore the equilibrium. This implies 
that rotational energy is converted into thermal energy and an 
enhanced neutrino emission is produced by a departure from the 
beta equilibrium. Thus, this mechanism can keep old millisecond 
pulsars (MSPs) warm, at temperatures ~ 10 5 K, which implies 



that the surface temperature of the MSP J0437-4715 should be 
20% below its measurement (Fernandez & Reisenegger 2005). 

On the other hand, cooling curves usually consider the ef- 
fects of nucleon superfluidity on the thermal evolution of NSs. 
Superfluidity is produced by Cooper pairing of baryons due to 
the attractive component of their interaction, and it is present 
only when the temperature T of the matter falls below a criti- 
cal temperature T c . The physics of these interactions is rather 
uncertain and very model-depende nt, and so is the critical tem - 
perature obtained from theory (see lLombardo & SchulzefeOOll) . 
An important microscopic effect is that the onset of superfluidity 
leads to the appearance of a gap A in the specttum of excitations 
around the Fermi surface. This considerably reduces the neu- 
trino reactions and the specific heat involving superfluid species 
(lYakovlev et al.l l200lh . and additionally opens new neutrino 
emiss ion processes, name ly pair breaking and formation reac- 
tions (Flo wers et al.ll976l) . Taking these effects into account, NS 
cooling has been used to constrain the amplitude of the energy 
gaps by comparing predictions with the s urface temperatures 
measu red from young neutron stars (see Yakovlev & Pethic 
l2004l for a complete discussion). In particular. lPage et al. I (I200 1 
considers the minimal cooling paradigm, where no direct Urea 
processes are allowed and the cooling is enhanced by Cooper- 
pair emission processes. They find this mechanism to be consis- 
tent with observations as long as the critical temperature 3 f2 of 
the neutrons covers a range of values between T™ n < 0.2 x 10 9 
K and T™ x > 0.5 x 10 9 K in the core of the star. 
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In this paper, we include the effects of superfluid- 
ity in rotochemical heat ing, buildin g on the framework o f 
Fernandez & Reiseneggei] d2005l) and lReisenegger et alj ([2006). 
The order of mag nitude of these effects was estimated in 
Reiseneggei] d 1997b . by considering the reactions to be totally 
forbidden, until the chemical imbalances r\ np ] = jj„ - fi p - fi/ 
(I = e,yu) exceed a combination of the energy gaps of either 
A„ + A p when direct Urea is allowed or A„ + A p + 2A< when 

only modified Urea operates, where A< = min (A„, A p ). This ef- 
fect lengthens the equilibrium timescale and raises the surface 
temperature. 

We consider how exactly the neutrino reactions are sup- 
pressed by super fluidity in the core, by avoiding the step-like 
approximation of lReiseneggerl (1 19971) . We calculate instead the 
reduction in the net modified Urea reaction rate and the emis- 
sivity in the regime where the energy gaps, the chemical im- 
balances, and the temperatur e are a ll relevant quantities. In this 
scenario, IVillain & Haen sel (2005) numerically computed the 
phase-space integrals of the net reaction rate for direct Urea and 
modified Urea processes, finding that Cooper pairing does not 
strongly inhibit these reactions when the energy gaps are not too 
large compared to both the temperature and the chemical imbal- 
ances. 

The structure of this paper is as follows. In Sect.|2]we present 
the basic equations of rotochemical heating and describe the su- 
perfluid input to compute these numerically. We explain our ap- 
proach to calculating the reduction factors and how they behave 
in the regime of interest. In Sect. [3] we describe our results and 
compare our prediction with the like ly thermal emission of PSR 
J0437-4715 dKargaltsev et al.ll2004l) to constrain the values of 
the energy gaps. We summarize our main conclusions in Sect.|4] 
Finally, in Appendices lAlandlBlwe explain in detail our numeri- 
cal and analytical approaches to computing the net reaction rate 
and emissivity. 



2. Theoretical framework 

2.1. Rotochemical heating: basic equations 

The framework of r otochemical heating used in this w ork is de- 
scribed in detail in IFernandez & Reisenegger (2005). Here we 
just point out the fundamental equations and the modifications 
required to introduce the effects of superfluidity. 

We consider the simplest model of a neutron star core, com- 
posed of neutrons, protons, electrons, and muons (npefi matter), 
ignoring the potential presence of exotic particles. 

The internal temperature, redshifted to a distant observer, 
T x , is taken to be uniform inside the star because we are model- 
ing the thermal evolutio n of a MSP over tim escales much longer 
than the diffusion time (Reisenegger 1995). Thus, the evolution 
of the internal temperature for an isotherma l interior is given by 
the thermal balance equation dThornd 19771) 



rBF.OO 



l;), 



(i) 



where C is the total heat capacity of the star, the total power 
released by the heating mechanism, the total power emitted 
as neutrinos due to Urea reactions, Ly F '°° the total power emit- 
ted as neutrinos due to Cooper pair-breaking and pair-formation 
(PB-PF) processes, and L™ the power released as thermal pho- 



tons. We briefly discuss the quantities C and L: 
and Sect. 12.51 respectively. 
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in Sect. UA 



The amount of energy released by each Urca-type reaction is 
rfnpi — /Ai - Mp ~ W ( I = e < AO- Thus, we write the total energy 
dissipation rate as 



I CO OO A y~ CO A f 1 " 

H ~ ^npe ^ npe + T lnpu m npp.i 



(2) 



where Ar„ p / = Y„^ p i - T p i^ n is the net reaction rate of the Urea 
reaction integrated over the core (indicated by the tilde) involv- 
ing the lepton /. 

The photon luminosity is calculated by assuming black-body 
radiation L~ = 4no-{R°°) 2 T A S0O , where R 00 and T s> „ are the ra- 
dius and the surface temperature of the star measured from 
an observer at infinity, respectively. To relate the internal and 
the surface temperat ures, the fully accreted envelope model of 
iPotekhin etail (1 19971) is used. 

The evolution of the redshifted chemical imbalances, which 
are also uniform throughout the core, are given by 



inpe 

■ 00 

^lnpp 



-7 AF 

t^npe 1 -* 1 npe 



■ 7 AF 

t^np 1 -* 1 npfi 



■ 2W ttpe Ml 



(3) 
(4) 



and W„ P p are constants 



where the terms Z„ p , Z npe , Z„ PI1 , W npe . 
that depend on the stellar structu re and are unchanged wi th re 
spect to their latest definition in Reisenegg er et al.l d2006t) , and 
Qfi is the product of the angular velocity and its time derivative 
(proportional to the spin-down power). 

The key new contribution in this work is the recalcula- 
tion of Ljf and Ar^ne, which diff er substantially from those of 
IFernandez & R eisenegger (2005), beacuse superfluidity strongly 
inhibits these reactions. 



2.2. Cooper pairing 

In the core, neutrons are believed to form Cooper pairs due to 
their interaction in the triplet 3 Pi state, while protons form sin- 
glet l So pairs. In addition, neutrons in the outermost core and 
inner crust are believe d to form singlet-state l Su pairs (type A) 
dYakovlev et all 1200 lb . The 3 P 2 (type B and C) state descrip- 
tion is rather uncertain in the sense that the energetically most 
probable state of ««-pairs (\mj\ = 0, 1,2) is not known, being 
extremely sensitive to the stil l unknown nn-interaction (see, e.g., 
lAmu ndsen & 0stgaard 1985). Taking this classification into ac- 
count, IVillain & Haensell (120051) solve numerically the suppres- 
sion caused by each type of superfluidity of the net reaction rate 
for direct Urea and modified Urea reactions out of beta equilib- 
rium, finding that the suppression caused by type A superfluidity 
is in between the suppression by anisotropic types \m j\ = (type 
B) and \ntj\ = 2 (type C) superfluidity, respectively. For sim- 
plicity, we consider the energy gaps for the neutrons A„ and the 
protons A p at zero temperature, redshifted to a distant observer, 
as parameters that are isotropic ( l So pairs) and uniform through- 
out the core of the NS. 

The phase transition for a nucleon species into a superfluid 
state takes place when its temperature falls below a critical 
value T c . This temperature is related to the energy gap at zero 
temperature A(T = 0); for the isotropic pairing channel 'So, 
A(r = 0) = \.16AkT c . In addition, when the transition occurs 
the amplitude of the energy gap depends on the temperature by 
means of the BCS equation ( Yako vlev et al.l20'0l1). which can be 
fitted by the practical formula of Levenfish & Yakovlevl fl994) 
for the isotropic gap given by 



5 = 



ACT) 
kT 



Vl -r/7^1.456- 



0.157 1.764 
+ 



^fjr c t/t c 



(5) 
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where 6 is the variable used in the phase-space integrals in 
Sect. 12. 3 1 that depends on the input parameters A„ and A ;; used 
throughout this paper. It is straightforward to check that the lim- 
iting cases are reproduced by the Eq. ([5]), i.e., 6 = when 
T _ T c and 6 _ A(T = )/kT when T <K T c . In addition, 
iLevenfish & Yakovlevl (1 19941) claim that intermediate values of 
T/T c are also reproduced by this formula with a maximum error 
smaller than 5%, which is accurate enough for the purpose of 
this work. 

Having defined the energy gap of the nucleon A, with i = 
n, p, it is possible to express the momentum dependence of the 
nucle on energy 5(Pj ) near the Fermi level, i.e. \pi - pp.\ <K pp., 
as in I Yakovlev et al.1 (l200lh as 



£i(Pi) = Pi ~ y]v 2 Fj (Pi- Pf;) 2 + if Pi < Pf„ 

diPi) = Pi + ^v 2 F ipi- p Fi ) 2 + A 2 if pi > p Fi , (6) 

where pp., vp., and p/ are the Fermi momentum, the Fermi veloc- 
ity, and the chemical potential of species i = n,p, respectively. 

2.3. Neutrino emissivity 

The most rapid reactions in NS cores are the direct Urea pro- 
cesses. However, as already mentioned, we assume in this work 
that these reactions are forbidden because direct Urea with su- 
perfluidity might drastically change the behavior of rotochemi- 
cal heating, leading to new conclusions that will be discussed in 
a forthcoming work. Additionally, these reactions can be kine- 
matically forbidden for a wide range of EOSs and central densi- 
ties. If this were the case, the so-calle d modified Urea reac tions 
(or Murca reactions) would prevail ( Yak ovlev et al.l [200 ll) be- 
cause they involve an additional spectator nucleon N that allows 
energy and momentum conservation. In general, these reactions 
are 



n+Nj^p + Nf+e + v e 
p + Nj + e~^>n+Nf + v e , 



(7) 
(8) 



where the subindices i and / represent the for its initial and final 
states of the spectator nucleon. If this nucleon is a proton, these 
reactions are called the proton branch of Murca, and if it is a 
neutron, they are called the neutron branch of Murca. 

We write the neutrino emissivity and the net reaction rate due 
to Murca reactions involving the lepton / and integrated over the 
core, respectively, as 



AT, 



npl 



J Jtl T~-8 . J jP T--8 

£ 'M,l - oo ■ £ l M,T °°> 



L nl T n 



in T 7 , P l tP t 1 

M,Y OO ' • l 7lJr* n 



(9) 
(10) 



where constants L„i and L p i are defined in terms of the neutrino 
luminosities for a nonsuperfluid NS in beta equilibrium, as 



i eq r 

l oo Jcore 



47Tr 2 e A S a (n)e-"' i 'dr, 



6* 



(ID 



where a indicates both the bra nch of the Murca process and th e 
lepton involved in the reaction (Fernan dez & Reiseneggerl2 005). 
The term S a is a slowly varying fun ction of the baryon number 
density n (e.g jYakovlev et ai]|2001|) . and A and (p are the usual 
Schwarzschild metric terms. The quantities 1^ and are di- 
mensionless phase-space integrals that contain the dependence 
of the emissivity and the net reaction rate, respectively, on the 
chemical imbalances rj™ pl and on the energy gaps A„ and A p . 



To introduce these integrals, it is useful to define the usual 
dimensionless variables normalized by the thermal energy IcT, 
of 



kT 



e v , . Vnpl 

Xy = — , and ti = , 

kT h kT 



(12) 



which represent the energy of the nonsuperfluid degenerated par- 
ticle j, the neutrino, and the chemical imbalance involving the 
lepton /, respectively, while for the superfluid nucleon i we write 



Xi = XF '^ P ' kT — and Z< = SgnOO sjxf + tf- 

In terms of these variables, 
60480 



(13) 



jN _ 
1 M,e — 



q y"»CO r\OQ r\00 r\00 r\OQ r\00 

— ■ J dxyxl J J J J J dx n dx Nj dx Nf 



115137 

dx p dx e f(z„ )f(ZN, )f(zN f )f(z P )f(x e ) 

X \d(x v +%i-Zn~ Z Ni - Z Nf -Zp- X e ) 
+ 6(Xy -ijl-Zn- ZN, ~ Z Nf ~Zp- X e )j 



(14) 



and 



60480 



Q y-»CO y-»CO y-»CO y-»CO r\00 /-*00 

— ■ J dxyxl J J J J J dx n dx N .dx A 



11513z 

dx p dx e f(z„)f(z Ni )f(zN f )f(Zp)f(x e ) 
x [d(x v z„ - Zn, ~ z Nf -z p - x e ) 
- 6(x v z n - z N , - z Nf - Zp- x e )\ , 



(15) 



where /(■) is the Fermi function f(x) = 1/(1 + e x ) and the nu- 
merical factor in front of the integral is to normalize it to 1 when 
the energy gaps and the chemical imbalances are zero. 

In the nonsuperfluid case (i.e., 5„ = 8 D = 0), these int egrals 



reduce to the polynomials calculated by Reisene gger! (1 19951) 



I N M,e(8n = Sp = 0) = FuiZi) = 

22020£ 2 5670£f 420^f 

1 + T + -r + 



9$ 



115137T 2 115137T 4 115137T 6 115137T 8 ' 

14680^, 7560^ 3 840£, 5 24^/ 



115137T 2 115137T 4 115137T 6 115137T 8 



(16) 



(17) 



which are the same for the neutron branch and the proton branch. 

Since the z-variables defined in Eq. ( TPI ) depend on the en- 
ergy gaps, the equality between I" Me and I p Me , or between I" MT 

and I 1 ', is no longer satisfied if the gaps are different. 



2.3.1. Reduction factors 

The phase-space integrals in Eqs. ([T4l and (TT3T > do not have an 
analytical expression when one or both of the reacting nucle- 
ons are sup erfluid. Thus, their calculation must be done numer- 
ically, as in I Villain & Haensel d2005l) . A natural way to account 
for the suppression produced by Cooper pairing is to define the 
so-called reduction factors as the ratio of these superfluid inte- 
grals to their non-superfluid limits 



(18) 
(19) 



In principle, to calculate these reduction factors, a five- 
dimensional integral needs to be computed, because only one 
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dimension can be eliminated by integrating out the electron vari- 
able in Eqs. (TBt and (15[ . On the other hand, if one of the nu- 
cleons is not superfiuid, one can eliminate more dimensions by 
integrating out the nonsuperfiuid variables. For instance, if we 
consider both the neutron branch and the protons as the only su- 
perfiuid species, we can integrate out the electron plus the three 
nonsuperfiuid neutrons, obtain a two- dimensional integral that 
has to be calculated numerically (see [Villain & Haensell 120051 
for the formulae and details). However, it may be a problem to 
calculate these integrals efficiently including the new superfiuid 
luminosities and net reaction rates in the rotochemical evolution 
equations. 

The main difficulties in performing the integration of Eqs. 
(fT4l i and (O are 

- infinite sizes of the integration domains; 



Villain & Haensell |2005) calculate these integrals numeri- 
cally via the so-called Gauss-Legendre quadrature, using loga- 
rithmic variables scaled to external parameters r\\, 5 n , 6 p . In this 
way they cover a wide range of the infinite integration domain, 
for the parameters in Eq. ( fT9b in the range 6n e [0, 10 3 ] and 

6 e [o, io 4 ]. 

We require a fast integration method because we need to 
evaluate eight integrals (emissivity and net reaction rate for two 
nucleon branches and two leptons branches) for each time-step 
of the thermal evolution. F or this reason, we u s ed a s lightly dif- 
ferent method than that of Villain & Haense I (120051) . although 
we also implemented their code to calibrate ours. We chose the 
Gauss-Laguerre method, which is accurate enough for a few 
evaluations when the integrands are asymptotically exponen- 
tially decaying functions, as is the case for the Fermi functions 
(see AppendixlAlfor a detailed explanation). 

One striking feature of the reduction factors is that, for rel- 
atively high values of the energy gaps and chemical imbalances 
compared to th e thermal scale, they tend to be independent of 
the temperature. Villain & Haensel (2005), considering only one 
superfiuid particle species, claim by graphical inspection that, 
when 77 > A and A > lOkT, the reduction factors become func- 
tions of A/77 only. By using our numerical results, we show in 
Appendix B that this is a good approximation if A > 30kT and 
77 > A. In the limit of zero temperature and in the presence of 
one superfiuid variable we find the analytical expressions (see 
AppendixlBlfor more details) 



28(A/77) 2 + 105(A/77) 4 + ^(A/77) 6 + ^(A/t/V 



x In 



A/77 



1+ vi - (A/77) 2 



x |l + ^(A/^ + ^(A/,» 



105 



V 1 - (A/?7) 2 
4 : 1873 

105 



80 



(A/77)^20) 



fl M ,r(A/77) = I21CA/J7)' + —(A/77) 4 + -jj-CA/i/)' 



x In 



A/77 



_ . V 1 - (A/77) 2 

1+ Vl -(A/77) 2 j V 

x (l + ^(A/77) 2 + ^(A/77) 4 + ^(A/77) 6 ) , (21) 



where A is the energy gap of one of the nucleons, i.e., the proton 
if we consider the neutron branch and the neutron for the pro- 
ton branch. It is straightforward to verify that these expressions 
satisfy the limiting cases HmA/ n ->aR = 1 and liiriA/n^i R — as 
expected. In addition, in the limit kT <§c 1 the functions in Eqs. 
( TToT l and ( fT71 ) tend to the highest power of the polynomial, i.e., 



~ 11513s 8 

expressed as 



and Hffl = 



24? 

115137T 8 



W.0 



9? 



115137T 8 

24£ 8 

115137T 8 



R e (6/g), and 



Thus, the integrals can be 

(22) 
(23) 



- external free parameters 77;, 6„, and 8 P , which can be very i.e 
large; 

- many integration dimensions (up to five). 



where, doing some algebra , one can verify that these approxima - 
tions satisfy the identity of Flores-Tulian & Reisenegger (2006), 
// f) _ 3/ r ^ i n AppendixlBl we compare these formu- 



lae with the numerical non-zero temperature calculations, find- 
ing that they agree to a very good approximation when A/T > 
30. 
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Fig. 1. Contour plot of the reduction factor in the net reaction 
rate in the neutron branch of Murca in the limit of zero tempera- 
ture. The numbers are log{R" Mr ). 



When two nucleons are superfiuid, we cannot find an ana- 
lytical expression, but we can simplify the problem to a three- 
dimensional, bounded integral. This method is obviously much 
faster than the finite temperature one with the quadrature inte- 
gration in Appendix [A] However, it is an approximation to the 
problem, and only works when the temperature is quite low com- 
pared to the chemical imbalance and energy gaps. Furthermore, 
if we analyze the case with two superfiuid nucleons, the limit 
works when the chemical imbalance overcomes a certain energy 
threshold A,/„- imposed by the energy gaps A„ and A p . If we con- 
sider the reaction n + n, — > rtf + p + e~ + v e at zero temper- 
ature, the neutron is energetically allowed to decay only when 
2/4, - 2A„ > fi n + A„ + fi p + A n + ji e , which imp lies the simple 
condition 77^ > 3A„ + A r (see lReiseneggerll997l for a schematic 
justification of this). Finally, the threshold imposed by the gaps 
A t hr for both branches of Murca reactions are 



AfAr 
A,Z, r 



3A„ + Ap 
A„ + 3A ;; 



for the neutron branch and, 
for the proton branch. 



(24) 
(25) 



In Fig.Q] we plot the reduction factor in the net reaction rate 
of the neutron branch of Murca under the approximation of zero 
temperature. As shown in this figure, there are no reactions when 
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77 < A t i,r — 3A„ + A p . Remarkably, the contour lines of the reduc- 
tion factor are almost straight lines that coincide with the slope 
of the contour levels of equation 3A„ + A p . This means that the 
reduction factor in this regime is close to a certain function of 
A t hr, which is valid for both branches and also for the reduction 
factor of the emissivity. 



0.6 

u 

DC 

0} 

DC 0.4 
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log(A„/T) = 3 
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1 1.2 1.4 1.6 1.8 2 

Fig. 2. Ratio of the reduction factor of the emissivity R e to the re- 
duction factor of the net reaction rate Rr for the proton branch of 
Murca reactions with superfiuid neutrons as a function of t]/A„ 
for different values of A n /T. Ratio of the zero temperature ap- 
proximations given by Eqs. ( l20b and (fZTJ is given by the black 
dash-dotted line. 



Since we obtain the reduction f actors in the net reactio n 
rate that are similar to those found by Villa in & Haensell (|2005), 
these can be seen in their work. The behavior of the reduction 
factor in the emissivity is also similar to that in the net reaction 
rate. Thus, to illustrate the difference between these two quanti- 
ties, we plot their ratio R^/Rr in Fig. [2] for the range where the 
chemical imbalance slightly exceeds the energy gaps for differ- 
ent values of the temperature compared to the gaps. We show in 
Sect. [3] that this is the regime of interest for rotochemical heat- 
ing. We can verify using the formulae in Eqs. d20"b and (f2TT> that 
the ratio R € /Rr approaches a value of one when 77/A — > 00. The 
ratio also decreases as the temperature diminishes compared to 
the energy gap until the limiting case of zero temperature. This 
is even more dramatic if we have two superfiuid particle species. 
In conclusion, for very low temperatures relative to the other en- 
ergy scales and chemical imbalances slightly above the energy 
gap, the suppression of the emissivity is much higher than that of 
the net reaction rate. Physically what happens is that the neutrino 
released in the reaction escapes with a very small amount of en- 
ergy, which at zero temperature is proportional to the excess of 
energy ~ (77 - A,/,,-) that can be arbitrarily small as 77 approaches 

Affcr. 

We can finally incorporate the recalculation of the emissiv- 
ity and net reaction rate into the rotochemical equations when 
one of the nucleon species is superfiuid, by applying the Gauss- 
Laguerre quadrature method to the entire evolution. In contrast, 
when two superfiuid nucleons are present, we separate the evo- 
lution into two regimes and follow their respective approaches. 
While 77 < A t hr, we use the Gauss-Laguerre quadrature method 
explained in AppendixlAl which happens for the early evolution 
of rotochemical heating (see Fig.0. When 77 > A,/ lr , we use the 
analytical approximation of the integral presented in Appendix 

E 



2.4. Specific heat suppression 

When T reaches T c , there is a discontinuous increase in the spe- 
cific heat that is characteri stic of a second order phase transi- 
tion. lYakovlev et al.l d2001l) claim that this increase is by a fac- 
tor of ~ 2.4 with respect to the normal-matter specific heat. 
Subsequently, when T <k T c , an exponential-like suppression 
occurs because of the gap in the energy spectrum. In practice, 
these effects are taken into account by using control functions 
that multiply the unpaired values of the specific heat at con- 
stant volume Cy. These have been c alculated by several au thors 
for various tem perature regim es, e.g. lPizzochero et all d2002l) for 
T <K T c and lMaxwelll dl979h for Q.2T C < T < T c . The former 
set of authors obtain an exponential decaying control function of 
the form -^e~^ kT and the latter, which is the result we adopt, 
obtain the function 



C v ,n x 3.15 



C ' N „-l.l6T c . N /T 



T I T 
2.5- 1.66 +3.64 

T C ,N \ T C .N 



(26) 



where T c ^ is the critical temperature of the nucleon N and Cyjv 
is the specific heat of the norm al matter case, as defined in 
iFernandez & Reiseneggerl (120051) . As might be noticed, t his ex- 
pression also reproduces the result of iPizzochero et al. (2002) 
in the low-temperature r egime. Th e value at the lower limit of 
the expression found by iMaxwelll d 19791) is indeed Cyt{T - 
Q.2T c )/Cv,n = 0.005, which is sufficiently small to ensure a the 
lepton contribution to the specific heat dominant. This formula 
also represents the discontinuous increase at T — T c , which 
implies that C\ "^{T = T c )/Cv,n = 2.42, in agreement with 
lYakovlev et al.l d2001l) . Finally, an important remark is that the 
minimum value that the specific heat can reach is the leptonic 
contribution, which is Cy = Cy, e + Cy^. 



2.5. Cooper pairing emission 

Another feature of the superfiuid state is the appearance of new 
neutrino reaction mechanisms due to the Cooper pairs. These 
are the Cooper pair breaking and pair formation processes pro- 
posed bv lFlowers et al.l (Tl976). These authors claim that a super- 
fluid neutron star can be considered as a two-component system, 
which consists of paired quasiparticles and elementary excita- 
tions above the condensate. Their associated quasi-equilibrium 
densities are controlled by the processes, which are prevalent at 
temperatures close to T c and successively suppressed at lower 
temperatures because of an exponential decrease in the number 
of unpaired particles. Schematically, these neutrino reactions are 



NN 
N + N 



N + N + v + v 
NN + v + v 



pair breaking (PB), (27) 
pair formation (PF), (28) 

where AW denotes the Cooper pair and N the excitation. 
These authors find an emissivity Qpb,pf ~ 10 28 (/c77[MeV]) 7 
[erg cm -3 s -1 ] when T < T c . Thus, these reactions could cer- 
tainly affect the thermal evolution of the early stage, since the 
order of magnitude of the emissivity for the Murca processes is 
Qmu ~ 10 29 - 30 (/cr/[MeV]) 8 [erg cirT 3 s -1 ]. However, the rele- 
vant question here is whether or not it affects the late stage, i.e., 
when the photon luminosity is the dominant cooling mechanism. 
In this sense, we estimate the PB-PF emissivity for a a typical 
scenario of rotochemical heating, where the temperature that it 
reaches in the late quasi-steady state is kT ~ 10~ 3 MeV (see 
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iFernandez & Reisen egger 2005). In this limit, in which T <*c T c 
dFlowers etal.ll 19761) . 



Qpb,pf~ 1028 (^7) VaA?V 2AAT erg cm- 3 s- 1 , (29) g 



which clearly becomes negligible when the values of the en- 
ergy gaps are relatively large compared to the thermal scale, say 
A/kT ~ 10 2 - 10 3 , because the exponential behavior rapidly 
suppresses the effects of these reactions. For A > 0.01 MeV, we 
check that the PB-PF processes become irrelevant compared to 
the photon luminosity, and therefore do not contribute to the total 
luminosity at this stage. 



3. Results and discussion 

3.1. Evolution 

We show the results of numerically solving the rotochemical 
evolution Eqs. (Q}, ®, and ©, considering the inputs in Sect. 

m 

Hereafter, we model the neutron s tar structure by the A18 
+ Sv + UIX* equation of state (EOS) dAkmal et al.lll998l) . The 
most relevant feature of this EOS in this work is that it allows 
direct Urea reactions for electrons above a density po = 1.59 x 
10 15 g crrT 3 , which corresponds to the central density of a 2 M 
star. On the other hand, the threshold for direct Urea reactions 
of muons lies at a higher density in a non-causal regime. We 
consider stars below this mass limit, in which no direct Urea 
processes occur. 

The results are shown in Fig. [3] where in the upper panel 
we plot the solution to this set of equations for the nonsuper- 
fluid case, and in the lower we compare these to our results con- 
sidering superfluidity of neutrons for a constant energy gap of 
A„ =0.1 MeV by integrating numerically the emissivity and net 
reaction rate without any approximation. 

Except for the superfluidity energy gap of the neutrons in the 
lower panel, both panels in Fig. [3] have the same input parame- 
ters. This figure indicates that during the neutrino cooling stage, 
the differences between the nonsuperfluid and superfluid cases 
are not noticeable. After 10 5 yr, the temperature drops more 
rapidly with time in the superfluid case because the specific heat 
of the neutrons is suppressed. The photon cooling, which starts 
to dominate in both cases, is unaltered by the presence of the 
Cooper pairs. 

Later on, the nonsuperfluid case reaches a quasi-steady state, 
in which the rate at which the spin-down modifies equilibrium 
concentrations is the same as the rate at which the reactions drive 
the system toward the equilibrium conc entration, with he ating 
and cooling also balancing each other (seelReisenegg erTl995l and 
Fernandez & Reisenegger 2005 for more details). The timescale 
on which the superfluid case reaches this quasi-steady state is 
longer because the Murca reactions are highly suppressed when 
the chemical imbalances are smaller than A„. But, immediately 
after the chemical imbalances overcome the value of the energy 
gap of the neutrons, several neutrino reactions are produced, 
which drastically reheat the star, causing a rapid increase in the 
temperature to finally reach the quasi-steady state. The subse- 
quent evolution can be approximated by the simultaneous solu- 
tion of Eqs. (Q}, ®, and © with their left-hand sides set equal 
zero, as discussed in Sect. 13.21 Finally, the chemical imbalances 
at this final evolution stage reach a higher value than in the non- 
superfluid case, lengthening its timescale to reach quasi-steady 
state. This implies that the star can store more chemical energy 
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Fig. 3. Evolution of the internal temperature T, 
temperature r s co and the chemical imbalances rj™ pe , rf^ pu 
star with the mass fi xed to the PSR J0437-4715, i.e. 1.76M 
dVerbiest et al.ll2008h . built with the A18 + 6v + UIX* EOS, 
with the initial condition = 10 9 K and null chemical im- 
balances. The spin-down is assumed to be due to the magnetic 
dipole radiation, with dipolar field strength B = 2.8 x 10 8 G and 
initial period of Pq = 1 ms. The error bar is the 90% confidence 
level for the surface tempera ture measured for the PSR J0437- 
4715 dKargaltsev et al. 2004) at the current spin-down parame- 
ters. Upper panel: nonsuperfluid case (null energy gaps). Lower 
panel: superfluidity of neutrons with A„ = 0.1 MeV (dashed 
line). 



and it is released it later in the time evolution compared to the 
nonsuperfluid case. Therefore, the heating when Cooper pairing 
is present is more efficient to ensure that the MSP at higher tem- 
peratures during the quasi-steady state, compared to the normal 
matter case. Moreover, the choice of superfluid of neutrons with 
A„ =0. 1 MeV and nonsuperfluid protons can explain the 90% 
confi dence level of the surfa ce temperature of the MSP J0437- 
4715 (Kargaltsev et al. 2004), which the nonsuperfluid case can- 
not. 

This aforementioned effect of su perfluidity on rotoc hemi- 
cal heating was already predicted by Reiseneggej] (1 19971) via a 
rough estimation; this author claimed that the neutrino reactions 
opened when they overcome the threshold A,/,,., that for Fig. [3] 
is A,/„- = A„ (see Sect. 12.3. Il l, ignoring the temperature depen- 
dence of these reactions. By looking at in Fig.|3]and solving the 
evolution equations for several combinations of A„ and A p in the 
range of A,/„ = 0.05 - 1 MeV, we verify that for the Murca p ro- 
cesses, the zero temperature approach of lReiseneggeri (1 19971) is 
valid, because in the quasi-steady state T] npe ,r] nptl > A,j lr . To a 
good approximation, this means that the net reactions rates and 
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the emissivity indeed do not depend on the internal temperature, 
as we showed in Sect. I2.3.TI 



3.2. Quasi-steady state 

Sufficiently old stars will have reached the quasi-steady state, 
where t x = 0, r\™ pe = 0, and f}™ pft = 0. In this state, the Eqs. (D, 
(0), an d © can De written as 



3.2.1. Analytical approximation 

The net reaction rates AT„ pe and AT npil are determined entirely 
by Eqs. (l3TT l and (|32] >. Solving these equations, we obtain 



AT npe = 2QQ Z " pf,W " pe Z "" W " P " = -K e Dh, 

" 7 7— 72 

^npe^npp ^ n p 

Af np „ = 2 ^ Z n P eW npll -Z np W npe 

,l Ff* 7 7 72 " 

pe^n pp p 



(34) 
(35) 



2W npe Qh 



npe 

7 Af 

*- , npe Ly * npe 



npfi 



2W OD = 7 AF 



- 7 Af 



(30) 
(31) 
(32) 



where the superscript ^5 stands for quasi-steady and we have ne- 
glected the neutrino contribution of the Cooper pairing emission 
as explained in Sect. 12.51 This system of equations totally spec- 
ifies the temperature and chemical imbalances for a star with a 
certain value of QO. Figure [4] shows the solution to this equa- 
tions when only the neutrons are superfluid, such that A t u r = A„. 
This illustrates that the chemical imbalance always stabilizes at 
a value slightly higher than A,/, r , 



where the constants K e and K p depend exclusively on the 
stellar structure, i.e., th e EOS and the stellar mass (see 
[Fernandez & Reisen egger 2005 for the definition of these con- 
stants). 

In the quasi-steady state, in which > 10 2 , the polyno- 
mials in Eqs. ([Tol l and ( fTTI i are in a very good approximation 

9f 8 24£ 7 

= iisik? and = 11 vi -La * respectively. Thus, from 



the definitions (0, ( fTOb . ( fT8l ). ( 1191 ), and doing some algebra we 
substitute the expressions of Eqs. (f34t and d35l > to Eq. (f30b to 
obtain 



00 ,qs 
tfnpe 



oo,qs 
' ^Inpp 



min ( A„ + 3A p , 3A„ + A p j . 



(33) 



j oo.qs 



oo,qs 



n°°' qs K 
'/npe 



1 - 



1 - 



9 ^ne n Me ylnpe ) u pe^M,e ylnpe ) 

24 J E>« ( °°'1 S \ _i_ r nP ( 

'-'tieUfflY ylnpe ) + ^pe K M,r ylnpe I, 

r, 1 P" (r?°' qS \ J. J T> P {r,°°> qS \ 



24 f P'l /'n°°'' ? 'A J- f P P /n 00 ' 9 ^ 
Lj nfl lx M,T ylnpp ) " t " ^PP n MX ylnpp ) ) 



QQ.I (36) 



= 0.1 [MeV] 



where we have just rearranged the quasi-steady Eqs. ( l30l l. (f3~Tb . 
and ( l32l l. and no approximation has yet been made. 




Fig. 4. Chemical imbalance of the electrons in the quasi-steady 
state with only neutrons as superfluid particles as a function of 
OQ for three values of the energy gap A„ (solid lines) and the 
nonsuperfluid case (dashed line) for a 1.76M Q star calculated 
with the A18 + 6v + UIX* EOS. The horizontal lines indicate 
the values of the energy gaps: 0.01,0.05, and 0.1 MeV (dotted 
lines). 



We can also observe from Fig. |4] that the solution becomes 
closer to A tnr as we increase its value. For higher values of the 
gap, this is because the chemical imbalances at the quasi-steady 
state are higher. For a given spin-down rate, the net reaction rate 
is then fixed and the nonsuperfluid reactions become more effi- 
cient, increasing as ~ if . Hence, the reduction factors need to 
block more reactions, which happens for values of the chemi- 
cal imbalance closer to A tnr . We show in Fig. IB. 1 1 ( Appendix IBl 
that the reduction factors indeed become smaller as the chemical 
imbalance approaches the energy gap from above. 



log I HO I [radVs 3 



Fig. 5. Surface temperature predicted for the quasi-steady state 
for several values of A,/, f = min^A„ + 3A p , 3A„ + A p ) (solid 
lines) and the nonsuperfluid case (dash-dotted line) as a func- 
tion of the spin-down rate QQ with the same stellar parameters 
as in Fig. [4] 



The first approximation that we make, as previously dis- 
cussed, is to consider that the chemical imbalances are rff p f s « 
A°? . This limit is valid when the energy gaps are rela- 



Vnpp 



tively large, as we illustrated in Fig. 3] We checked this assump- 
tion for several solutions of the rotochemical heating evolution 
equations and found that it is a reasonable approximation for 
A thr > 0.05 MeV. 

Our second approximation is to neglect the ratios of the re- 
duction factors /Rm T . As we discussed in Sect. I2.3.T1 the 
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Fig. 6. Values of K e (solid lines) and (dash-dotted lines) for 
different EOSs for a range of masses from 1 M to the mass at 
which direct Urea with electrons is opened for each EOS. 



emissivity is more suppressed by the gaps than the net reaction 
rate when r\ > A. It becomes an acceptable limit when the values 
of the gaps are relatively large compared to the thermal energy, 
say A/kT > 10 2 , which is the limit of interest for rotochemical 
heating. We checked that when we consider A r /, r ~ 0.05 - 1 MeV, 
the ratio is R N M JR N Mr < 0.1. 

The first simplification tends to increase the predicted lu- 
minosity, while the second approximation tends to decrease it. 
Thus, both effects together tend to balance each other. Finally, 
placing both approximations together, we obtain the simple ex- 
pression for the bolometric luminosity of 



(K e + K^)A lh ,pQ.\ 



(37) 



whose physical interpretation is that the spin-down compression 
sets the number of reactions per unit time and each one of these 
reactions releases an a mount of energy A t i ir to reheat the star, 
as in the description of lReiseneggerl dl997l) . The approximation 
is even more accurate for the surface temperature, since 7^ oc 

(00 Cis\ 
L y j and the relative error in the lumimosity will lead to a 

smaller relative error in the surface temperature. 

In Fig. [5] we show the results of the predicted surface tem- 
perature in the quasi-steady state as a function of the spin-down 
rate by solving the equilibrium equations without approxima- 
tions and using the black-body law to relate L^' qs and T^, 
From this plot, we can infer that the curves are parallel and the 
surface temperature depends on some power of the spin-down 
that differs from the power of the nonsuperfluid case, which 
is higher. We can explain this behavior using the approximate 
expression oc |QQ| 1/4 , while for the nonsuperfluid case, 
iFernandez & Reiseneggerl d2005l) obtain rf M cc |Qii| 2/7 . 

In Fig. [6] we plot K e and K M as function of the radius for 
two s ets of EOSs. First, we show the two more realistic EOSs 
from lAkmal et al.1 (1 19981) for the core (A 18 + 6v and A18 + 
6u + UIX*), supplemented w ith that of iPethick et al.1 d 1995b 
and lHaensel & Pichonl (1994) for the inner and outer crust, 
respectively. Second , we plot four representative EOSs from 
Prak ash et al.l d 1988b . which open direct Urea reactions for stel- 
lar masses higher than 1 M . For this set of EOSs, we show in 
Fig.|6]that K e + a (2 - 6) x 10 47 [s 2 ], which, using Eq. Q7), 



infers a bolometric luminosity of 



L ^ (1 - 4)Xl ° \MeVj 



(P- 



20 



P 3 



erg s 



(38) 



where P_20 is the period derivative in units of 10~ 20 and Puis 
is the period in milliseconds. Finally, we can express the surface 
temperature using the expression for the luminosity given by Eq. 
03 as 



K e + K u 



\47TO-(R°°¥ 



1/4 



1/4 



(39) 



where <x is the Stefan-Boltzmann constant. For the set of EOSs 
that we use, we obtain 



Tf x (5.7 - 6.6) x 10 : 



/ A,/, r \ 
\MeV/ 



!/4 ( p \ 



1/4 



K. 



(40) 



3.2.2. Constraints on the energy gaps 

To explain the thermal emission of PSR J0437, as we showed in 
Fig. [3j it is necessary to invoke superfluidity, and the required 
values of the energy gaps then fall in a theoretically interesting 
range. In Fig. [7] we c ompare the observation al allowed range of 
surface temperatures ( Kargalt sev et al.l2 004) with the theoretical 
predictions for different values of A,/,,-, using one particular EOS. 
In this case, we obtain the allowed range 

0.05[MeV] < min (A„ + 3A p , 3A„ + A p ) < 0.2[MeV]. (41) 

By adding two more EOSs (BPAL 21, BPAL 9) that forbid di- 
rect Urea reactions in the allowed mass range, we expand the 
constraint to 

0.05[MeV] < min (A„ + 3A p , 3A„ + A,,) < 0.45[MeV]. (42) 

These results depend of course on the assumption of spatially 
uniform and isotropic energy gaps. In principle, these are unre- 
alistic assumptions because the energy gaps depend on the local 
density and therefore on the radius of the star, and the neutrons 
in the core are expected to form anisotropic Cooper pairs (see 
Sect.r 



3.2.3. Condition for arrival at the quasi-steady state 

To estimate the time taken for the star to arrive at the quasi- 
steady state, we consider that the chemical imbalances grow by 
the effect of Q.Q. in a unimpeded way, because the reactions are 
highly suppressed, until they overcome the threshold imposed 
by the gaps. Thus, the chemical imbalance associated with the 
lepton / evolves as 



(43) 



until it exceeds the threshold imposed by gaps, at the moment 
the quasi-steady state is reached. Thus, we set the condition for 
arrival at the quasi-steady state as rQ pl = A,i, r , which can be seen 
to be a good approximation from the lower panel of Fig. [3] Then, 
integrating the expression d43l over time from the initial spin 
period f, to the current value of the spin period P and using the 
previous condition, we obtain the upper limit to the initial spin 
period of 



1 

P 1 + 4tt 2 \W, 



npl\ ) 



-1/2 



1 

pT 



1 

+ - 

4 



Hhr 



Mev) 



-1/2 



ms, (44) 
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Fig. 7. Surface temperature predicted in the quasi-steady state as 
a function of the the stellar radius for a star built with the A18 
+ 6v + UIX* EOS for several values of A t h r in units of MeV 
(solid lines) and the current spin-down rate of the MSP J0437. 
Vertical dotted lines indi cates the mass range 1.76 + 0.2 M 
allowed for MSP J0437 dVerbiest et all [2008b . Dashed-dotted 
lines are the 90% conf idence contours of the black-body fit of 
Kargaltse v et alj (12004) to probable thermal emissi on from this 
pulsa r corrected for the latest distance of 157 pc dDeller et alJ 
2008). 



where P ms is the current value of the spin period in milliseconds. 

For the case of PSR J0437, in particular, we conclude that the 
initial spin period to reach the quasi-steady state with the upper 
limit of the constraint in Eq. (l42l . i.e., A t h r = 0.45 MeV, should 
be Pi < 2.7 ms. This value could in principle be compared with 
the initial period constraints from the cooling age of its white 
dwarf companion. Howev er, this constraint is highly uncertain 
dHansen & Phinnevll 19981) . and we are therefore unable to draw 
definite conclusions. 
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4. Conclusions 

We have studied the rotochemical heating of millisecond pulsars 
with only modified Urea reactions in the presence of uniform and 
isotropic Cooper pairing gaps for one or two nucleon species. 

Based on these assumptions, we have found that the chem- 
ical imbalances in the star grow to the threshold value A,/, r = 
min (A„ + 3A p , 3A„ + A ; ,), which is higher than in the quasi- 
steady state achieved in the absence of superfluidity. Therefore, 
the superfluid MSPs will take longer to reach the quasi-steady 
state than their nonsuperfluid counterparts, and they will have a 
higher luminosity in this state, given by 



C^(l-4)xW^) 



-20 



P 3 



ergs" 1 . (45) 



The constraint that we found for the energy gaps using our 
predicted effect i ve tem peratures and the black-body fit of 
iKargaltsev et"ai] (I2004T) to the UV emission of PSR J0437-4715 
is 

0.05[MeV] < A thr < 0.45 [MeV]. (46) 

In this sense, rotochemical heating presents an interesting tool 
for constraining the superfluidity parameters. In a future paper, 
we will include the density dependence of the energy gaps via 
the predictions of different theoretical models. 
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Appendix A: Gauss-Laguerre integration method 

We describe the numerical method used in this work to calculate the pha se-space integrals i nvolve d in the net reaction rate and 
emissivity of the beta processes considering superfiuid nucleon species. Villain & Haensel (2005) showed that Cooper pairing 
reduces these reactions and that there is no analytical expression to compute this suppression. These authors, therefore, solved the 
phase-space integrals numerically using the so-called Gauss-Legendre quadrature, logarithmic variables, and cut-off values to cover 
a wide range of the unbounded integration domain. In this w ork, it was more conveni ent for us to use the method presented here 
because it involves fewer evaluations than the method used in Villai n & Haensel |2005) for a reasonable precision. 

The method we show is motivated by the behavior of the integrand for the beta processes, which decays exponentially with 
each variable because of the Fermi functions. This happens asymptotically even when some particles are superfiuid. For this reason, 
among the Gaussian quadratures, a natural candidate for an interpolation polynomial are the so-called Laguerre polynomials because 
these are defined on the basis of the continuous functions in [0, oo], whose inner product is defined with weight function W(x) = e~ x 
(lAbramowitz & Stegunll972h . If the function to be integrated is the product of an nth-order polynomial and an exponential function, 
this method is exact using an nth-order Laguerre polynomial in the Gaussian quadrature. Therefore, a necessary condition for 
this method to work properly for an integral JT f(x)dx is that f(x) ■ W(x) has to be a smooth function, such as a polynomial 
dDavis & Rabinowit3ll975l) . 

We present how this method has been applied to the phase-space integrals involved in the net reaction rate for the direct Urea 
process to save notation. An extension to the Modified Urea processes is straightforward by adding two more superfiuid variables 
to the subsequent expressions. Considering the limits in the positive domain of each integration variable the integral is given by 

XOO y-»CO y-»00 j 1 

dxyxl I I dx„dx p xi^ ^ fUnZ n )f(jpZp)x [f(xy + ^-j n Zn-jpZp)-f{x v -^-j n Z n -j p Zp)]> , (A.l) 
J ° J ° [j„=±lj,,=±l J 

where zi = ^jx 2 + 6 2 (i = n,p) and /(■) is the Fermi function. In what follows, we argue that this integral satisfies the Laguerre 

quadrature condition explained above, and we then express the numerical formula to compute it. 

The integrand is an exponentially decaying function of the neutrino variable, and thus, it satisfies the quadrature condition. The 
integrand for the nucleons has the shape of a multiplication of two Fermi functions. For instance, if we consider the neutron variable 
we have f(± \Jx% + 5%)-f(a+ yx%~+ <$f), where a depends on the neutrino variable, the proton variable, and the chemical imbalance. 
It is easy to see that the asymptotic behavior of this function is exponential by separating both cases and imposing the limiting case 
x„ » 6„, which is the relevant scale for this variable. Thus, in this limit /(- yjx 2 + 8 2 ) ~ e~ x " and /( ^jx 2 + 6 2 ) ~ 1, while the other 

Fermi functions have the opposite behaviors f(a - sjx% + 6%) ~ 1 and f(a + ^x 2 + 5%) ~ e~ x \ respectively. Therefore, the nucleons 
also satisfy the quadrature condition, which may, however, also be applicable to higher values of these variables depending on the 
values that a take. In this sense, higher values of a, which depends on the inputs 6„, 8„, and will make the quadrature less accurate 
for a given amount of point evaluations. 

An additional step in applying this quadrature method is to define a relevant scale for each integration variable, as the scale of 
the argument of the Fermi function / (x v ± £ - j„z n - jpZp) in the integrand of Eq. ( IA.U . as in Villai n & Haensel (120051) . For the 

neutrino variable, this scale is S Xr = %+6 n + d p , and for the superfiuid nucleons they are S Xn = y]2S„ + 5 p and S x = ^26 p +6„. 
The numerical formula is, finally 

n v n n 

lDT = XXX X X S ^^ S ^-f0^k)f(jp^)^[f(^+^-i^-jp^-f(^-^-i^-ip^% (A.2) 

;,=i ;„=i i„=i j„=±\ j„=±i 



where z, = Ax 1 + 6l and z, = x 2 + 6] l. The values x-, , x, , and Xi are the roots of the Laguerre polynomials of order n y , n„, and 

n p , multiplied by S Xy , S Xn , and S Xp , respectively . The values W,„, W,„, and Wi are the weight factors of the Laguerre polynomials, 
which can easily be obtained from tables or using their usual recursive formulae. 

In practice, this method is used for Modified Urea processes where there are four superfiuid nucleon variables and the integration 
dimensions increase to five. 

In Fig. IA. H we provide an example where the Laguerre quadrature is more accurate than the Legendre quadrature over a few 
evaluations. However, to achieve high precision the Legendre-like integration method is preferable because it can be used to perform 
many evaluations raising the number of point evaluations, while it is more difficult to compute the roots and weights of the Laguerre 
method using recursive methods to high orders of these polynomials. 



Appendix B: Analytical approximation 

In the context of rotochemical heating, most of the MSP's lifetime is in a regime where the chemical imbalance is much higher than 
the temperature = ^ > 100). On invoking superfluidity, the chemical imbalances increase until they become of the order of the 
superfluidity energy gaps. Thus, a relevant case to analyze and calculate in this work is the limit 



€ ~ &n,p = -rdr » I- 



(B.l) 
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log(T|/A) 



Fig. A.l. Relative error in the integral for the Legendre quadrature and Laguerre quadrature as a function of the number of evalua- 
tions considered in each dimension for the integral of the direct Urea net reaction rate with £ = 5 n — 6 P = 10. 



The first analysis is given for the simplest case: the direct Urea process. The idea is to extend this analysis to the most general and 
complicated case: modified Urea with two superfiuid nucleons. After integrating over the electron variable, the phase-space integral 
for the net reaction rate is 

dx Y x 2 v I I dx n dx p f(z„)f(Zp) X {f(x v + £, - z„ - z p ) - f(x Y - £ - z„ - z p )}, (B.2) 

J —CO tJ — CO 

where Zi = sgn(x,) Jx 2 + S 2 , with i = n, p. 

In the regime of interest, the Fermi functions behave almost lik e a step function, since t he thermal scale is negligible compared 
to the other relevant scales. Moreover, one of the conclusions of Villai n & Haensell (|2005), obtained by means of their numerical 
results, is that for large enough energy gaps, say log(A/r) > 1, and as soon as rj > A, the reduction factors of the net reaction rates 
for direct Urea and modified Urea reactions become very nearly independent of the temperature. This motivates us to take the limit 
of zero temperature. Thus, replacing the Fermi functions by step functions, i.e., f(x) — > ®(-x), the integral gives 

I D ,T = I dXyxl I I dx n dx p ®(-Zn)®(-Zp) X {®(-X Y - £ + Z n + Zp) ~ &(~X V + £ + Z n + Zp)}. (B.3) 

U0 *J — CO yj —OO 

Rearranging the limits of integration and eliminating the sgn(-) function from z, = sgn(x,) yjx^ + 6 2 i7 the integral can be reduced to 

Ids = f f Q f Q dx v dx ndx p x 2 v ® - x v - ^jx 2 n + 5\ - ^Jx 2 + 5*j . (B.4) 

By applying these approximations, the reactions can be seen to be open only when 5„ +6 P < (see lReiseneggerl 1 9971 for a schematic 
justification). 

The integral in Eq. ( lB.4t appears not to have a general, closed expression in terms of the three parameters ^,6„,6 P , but it is 
possible to perform two out of three variable integrations. By considering only one superfiuid nucleon, an analytical expression can 
be found. In what follows, we explain these results. 

The first step is to define an appropriate change of variables 



x 2 + 6 2 



u v = — and ui = -J — ~ 2 — , (i = n,p), (B.5) 

where the Jacobian of this transformation is J(u v ,u n ,u p ) = & , considering that A,,/j] = and the new 

V«S-(A»/'7) 2 yul-( A p/D 2 

limits of integration must satisfy the upper limit imposed by the step function and the lower limit defined by the change of variables. 
Thus, the new integral is 



r' r'-«» r l -""-"» , u„u p 

/or = £ I I I du n du p du v u v — — . (B.6) 

Ja„/„Jv„ Jo ^Jul-iAJtf^juj-iAp/tf 



Integrating over the neutrino variable, we obtain 



£ 5 f 1 f 1 """ , U n U p 

I D ,r = — du n dUp(l-u n -u p ) . (B.7) 
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The integration over the proton variable yields 

e5 /-1-V7 



ll„ 



In 



V«l - (AJn) 



:K r (u n ,A p /j]) 



(B.8) 



where 



K T (u n ,A p /rj) = - 



3(A /7 /7 7 ) 2 



(A p /i]f 



•In 



1 - Un + V(l " ««) 2 - (Ap/l]) 2 



Ap/r] 



(B.9) 



■^(1 - ««) 2 - (Ap/r/) 2 



(1 - u n ) + 



The latter integral in Eq. jB.9\ appears not to have an analytical closed formula. For this reason, the solution hereafter is taken 
numerically. The method for the emissivity integral is completely analogous, but with the only difference that we replace x 2 with x\ 
in Eq. ( IB. 4b . The formula obtained in this case is 



t6 ri-a, 

lD * = 4 

4 J&n/TI 



'l-A p /i7 



du n 



V«5 - (A„/r?) ; 



:K € (u n ,A p lrj) , 



where 



KMnAp/jj) = KAplijf 



-{A p lrff(\ -k„) + (1 -w„) : 



V(l-«„) 2 -(A p //7) 2 



In 

(Ap/7?) 2 



1 - K„ + -y/(l - M„) 2 - (Ap/l]) 2 



A P h 



(B.10) 



(B.ll) 



(1 - m„) 4 + [83(1 - u„f + 16(A P /^) 2 ] . 



The final step is to extend the same reasoning to the modified Urea processes, and is natural on the basis of the previous analysis. 
Here one has to distinguish between the neutron branch (superscript n) and the proton branch (superscript p), but for the similarity 
relation I P MY (6„,6 P ) = I n MT {S p , 6„) only one expression is needed. Thus, without loss of generality, we only write the neutron branch 
expression of the net reaction rate of Murca 



4r = — I I ( du„du n .du n 



yju\ - (AJtj) 2 - (AJii) 2 ^u 2 „ f - (AJij)- 



■K r (u n + u n> + u nf ,A p lrf), 



(B.12) 



where Kr(-, •) is the function given by Eq. ( |B.9t and n,- and rif represent the initial and final neutron spectator, respectively. In the 
same way, the emissivity has an expression in terms of the function K e (-, ■) defined by Eq. (IB. 1 U : 



l n Me = — I I I du„du ni du„ 

4 Ja„/ij JAJij JAJij 



V« 2 - (A„//,) 2 - (AJtj) 2 ^u 2 nf - (A„/t,) 2 



K e (u A P h). 



(B.13) 



In both cases, the net reaction rate and the emissivity, the condition to recover the non-superfluid case asymptotically is satisfied. 

The advantage of this treatment is that the numerical integration of the net reaction rate and the emissivity for the Urea processes 
is much faster than those without these approximations. First, the region of integration is bounded, while those of the initial integrals 
are not. This reduces considerably the number of points used in each integral. Secondly, the dimensions to integrate are reduced to 
one in t he case of direct Urea an d to three in the case of modified Urea. In practice, to compute these integrals we use the tanh-sinh 
method (Taka hasi & Morill 19731) . which integrates the singularities in A,Jt] and/or in A p /tj quite rapidly. 

The integrals in Eqs. ( IB. 81 ), ( IB. 10b . ( IB.12I I. and ( IB . 1 3b can be solved analytically when only one superfluid reactant is present, 
which is characterized by the energy gap A. The reduction factor R (for all the previous cases) obtained from this simplification is 
given by the integral 



R = (a + 1) I du- 

Jam 



u(\-u) a 

V« 2 - (A/7,) 2 ' 



(B.14) 



where a — 4 and a — 5 for Durca processes (net-reaction rate and emissivity, respectively), and a — 6 and a — 7 for Murca 
processes (net-reaction rate and emissivity, respectively). The solution to this integral is 



R(AJi 1 ) = P(A/Tj)ln 



1 



1+ Vl-(A//?) 2 J 



Q(A/r 1 )^\-(Alr 1 ) 2 , 



(B.15) 
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where P(-) and Q{-) are the polynomials 

P(x) = 10x 2 + ^x 4 and Q(x) - 1 + ^x 2 + ^x 4 for the Durca net reaction rate, (B.16) 
2 6 3 

45 15 97 113 

P(x) = 15X 2 + —x 4 + —x 6 and Q(x) = 1 + — x 2 + -3-x 4 for the Durca emissivity, (B.17) 
2 8 4 8 

P(x) = 21x 2 + ^^-x 4 + ^^-x 6 and Q(x) = 1 + -^x 2 + ^r— ' x 4 + — x 6 for the Murca net reaction rate, (B.18) 
w 2 8 20 40 5 

, 4 105 . 35 s 551 , 4327 , 1873 fi 

P(x) = 28x 2 + 105x 4 + —x 6 + —x 8 and Q(x) = 1 + — x 2 + ~^-x + ~oq~ x for the Murca emissivity. (B.19) 



- Laguerre quadrature 
■ Legendre quadrature 




Fig. B.l. Reduction factor R for the net reaction rate of direct Urea with one superfiuid particle as a function of the ratio 77/ A for 
several values of the ratio A/T (dashed lines) using the method detailed in Appendix|A] and the limit case 77 A — > calculated from 
Eqs. dB~T5l l and SbIM (solid line). 



Finally, we show in Fig. IB. 1 1 a comparison between our zero temperature approximation and the exact calculations for finite 
values of temperature, which demonstrates the accuracy of our approximation even for values of the energy gap compared that are 
not very large relative to the temperature, say A/T > 30. 



